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Abstract 

A  3D  simulation  tool  for  solid  oxide  fuel  cells  is  presented.  The  aim  of  this  work  is  to  predict  current  density,  flow,  temperature  and 
concentration  fields  in  order  to  compare  and  optimize  repeat  element  geometry  for  a  whole  stack.  A  commercial  CFD  tool  was  used,  solving 
mass,  momentum  and  energy  equations;  whereas  chemical  kinetic  equations  are  computed  from  external  sub-routines.  A  steady-state  case 
is  presented,  fed  with  hydrogen.  The  flow  is  laminar  for  both  air  and  fuel.  Radiative  heat  transfer  is  taken  into  account  between  inner 
surfaces.  On  boundaries,  convective  and  radiative  heat  transfers  are  assumed  at  external  surfaces  between  repeat  element  and  oven.  Due 
to  the  large  range  of  dimensions  (cells:  300  fxm  thick,  gas  channels:  1  mm  height,  whole  cell:  80  mm  x  80  mm)  a  fine  mesh  was  needed. 
Data  for  conductivities  and  kinetics  were  estimated  from  experiments  performed  in-house.  Simulation  results  are  presented  and  compared 
to  real  repeat  element  test  measurements  for  the  current-potential  characteristics. 
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1.  Introduction 

At  this  stage  of  SOFC  stack  development,  optimization 
of  geometry  is  almost  as  crucial  as  material  development. 
The  geometry  of  a  stack  is  complex,  with  manifold  and  gas 
channels  distributing  fuel  and  air  among  and  within  cells. 
The  solid  oxide  fuel  cell  implies  high  operating  temper¬ 
ature  and  introduces  issues  in  thermal  management.  For 
example,  the  thermal  behavior  will  influence  the  material 
stability  and  the  response  of  the  stack  to  fluctuating  inlet 
gas  compositions  is  not  commonplace.  As  experiments  are 
expensive  and  time  consuming,  there  is  a  need  for  tools,  to 
predict  current  densities,  temperature  field,  pressure  drop, 
etc. 

In  this  paper,  a  computational  fluid  dynamics  (CFD)  tool 
is  presented,  allowing  detailed  calculations  of  physical  phe¬ 
nomena.  The  accuracy  of  FLUENT™ commercial  software 
has  been  widely  demonstrated,  but  as  we  added  sub-routines 
to  calculate  reaction  rates  and  current  densities,  we  verified 
the  conservation  equations. 

Generally,  development  of  a  new  cell  structure  or  elec¬ 
trode  composition  is  performed  on  small  active  areas 
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(1  cm2).  For  the  stack  development,  larger  cells  are  needed 
and  power  densities  observed  on  these  (^50  cm2)  are  far 
less  than  those  on  small  cells.  We  decided  to  determine  cell 
resistance  and  polarization  losses  from  small  cell  measure¬ 
ments  and  compute  the  behavior  for  the  larger  geometry. 
Finally,  we  compared  the  simulation  results  with  experiment 
on  the  same  geometry. 

2.  Experimental 

Cells  are  anode-supported,  200  fxm  thick  with  a  6-8  \im 
thick  electrolyte  layer.  While  developing  cells  for  stack  test¬ 
ing  we  produce  mainly  two  sizes  of  cells.  On  the  one  hand, 
cells  composed  of  16  cm2  of  anode- supported  electrolyte 
(ASE)  with  1  cm2  active  cathode  area  are  performed  to  mea¬ 
sure  electrochemical  performances  of  new  compositions, 
and  for  quality  control  of  produced  cells.  On  the  other  hand, 
67  cm2  of  ASE  cells  with  50  cm2  active  area  are  produced  to 
build  stacks.  Both  sizes  are  produced  from  the  same  batches 
to  ensure  as  much  as  possible  identical  electrochemical  prop¬ 
erties.  In  the  present  study  we  fed  cells  with  hydrogen  as  fuel 
and  air  as  oxidant  to  compare  with  the  simulation,  where 
reforming  reactions  have  not  been  implemented  yet.  I  —  V 
characteristics  are  recorded,  at  constant  air  and  fuel  flux,  and 
for  several  temperatures. 
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Nomenclature 

Cij 

inertial  losses  coefficients 

Cp 

heat  capacity  of  the  mixture 

Dij 

diffusion  coefficient  for  species 
in  the  mixture 

Dj  ,i 

thermal  diffusion  coefficient 

Day 

Darcy  viscous  losses  coefficients 

Pj 

binary  mass  diffusion  coefficient 

Ef 

total  fluid  energy 

Es 

total  solid  energy 

F 

Faraday  constant 

F 

momentum  source 

AG0 

free  Gibbs  energy  of  formation 

hi 

enthalpy  of  i  species 

I 

unit  tensor 

j 

current  density 

Jleak 

leakage  current  density 

J, 

diffusion  flux  of  species  i 

kcff 

effective  thermal  conductivity 

Mout,  Min 

mass  flux  in  and  out 

Pi 

partial  pressure  of  species  i 

91 

heat  flux  through  boundaries 

n 

net  rate  of  production  of  species  i 

R 

universal  gas  constant 

^ohm 

ohmic  resistance  of  the  electrolyte 

tfpol 

resistance  for  polarization  losses 

sh 

heat  source 

T 

temperature 

Tout 

temperature  of  the  outlet  gas 

TJ cell 

potential  of  the  cell 

UNernst 

Nernst  potential 

TSocn 

cell  potential  in  open  circuit 

V 

velocity 

Unag 

velocity  magnitude 

Vi 

diffusion  velocity 

Xi 

mole  fraction 

Yi 

mass  fraction  of  species 

Greek  letters 

£ 

error 

y 

porosity  of  the  medium 

jl 

molecular  viscosity 

p 

density  of  the  mixture 

X 

stress  tensor 

V 

gradient 

2.1.  Small  cells  measurements 

The  aim  is  to  characterize  the  anode  supported  electrolyte 
(ASE)  to  define  the  model  parameters.  The  small  cell  is 
tested  in  near  ideal  condition,  with  massive  fuel  and  air  flux 
(fuel  utilization  less  than  10%),  and  platinum  and  nickel 
mesh  as  current  collectors.  I  —  V  measurements  have  been 
performed  from  690  to  890  °C.  Results  are  shown  in  Fig.  1. 


Fig.  1.  Characteristic  I  —  V  curves  for  1cm2  active  area  button  cells 
(16  cm2  ASE). 


H2 


Fig.  2.  Repeat  element  concept  with  internal  manifolding. 

2.2.  Single  repeat  element  tests 

Experiments  on  repeat  elements  are  based  on  the 
HTCeramix-LENI  concept  (see  Fig.  2),  with  internal  man¬ 
ifolding,  on  80  mm  x  80  mm  anode  supported  cells  with 
about  50  cm2  of  active  cathode  area,  with  gas  diffusion  lay¬ 
ers  on  both  electrodes  [1].  Interconnects  are  formed  with 
0.75  mm  thick  metallic  sheets  on  both  sides. 

The  present  paper  discusses  a  single  repeat  element,  so 
results  with  only  one  element  in  a  setup  are  given  here. 
The  setup  is  placed  in  an  electrically  heated  oven  with  a 
controlled  temperature  of  750  °C. 


3.  Theory 

Understanding  and  predicting  SOFC  behavior  implies  cal¬ 
culating  at  the  same  time  various  coupled  phenomena: 

(1)  thermal  behavior,  with  conductive,  convective  and  ra¬ 
diative  heat  transfer; 

(2)  multicomponent  gas  flow; 

(3)  electrochemical  reactions,  involving  partial  pressures, 
local  temperature,  ionic  and  electronic  conduction. 

Facing  those  issues  we  decided  to  use  the  commercial  CFD 
tool  FFUENT,  able  to  compute  the  two  first  points  above, 
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and  added  sub-routines  to  calculate  the  electrochemical  part. 
Electrochemical  behavior  of  the  cell  is  experimentally  fitted 
on  small  cells  where  we  could  partially  uncouple  the  phe¬ 
nomena,  and  then  applied  to  complex  geometry. 

3.1.  Assumptions 


porous  medium: 

d  _  _  - 

—  (pv)  +  V(pvv)  =  -Vp-\-  V(r)  +  F  (4) 

ot 

with 

f  =  /x[(V5  -1-  VuT)  —  f  Vu/]  (5) 


For  thin  anode  supported  cells,  we  assume  that  the  cur¬ 
rent  paths  are  straight,  perpendicular  to  the  electrolyte  layer 
[2].  At  high  temperature,  gas  mixtures  can  be  accurately 
represented  as  mixtures  of  ideal  gases,  so  heat  capacities 
are  functions  of  temperature  only  [3].  The  heat  capacity  of 
the  mixture  is  calculated  with  a  volume-weighted  mixing 
law  of  individual  heat  capacities.  Thermal  conductivities 
and  viscosity  are  calculated  from  a  mass- weighted  mixing 
law  of  individual  values.  Individual  values  are  calculated 
from  a  temperature- dependent  polynomial  curve.  The  cell 
thermal  conductivity  is  calculated  from  a  Ni-YSZ  anode 
cermet  formula  [4]. 

According  to  Reynolds  numbers,  between  1  for  fuel  flow 
and  10  for  air  flow,  the  flow  is  assumed  to  be  laminar.  Anode 
support,  anode  active  layer,  electrolyte  and  cathode  active 
layer  are  modeled  as  one  discrete  unit.  Reactions  are  defined 
on  each  surface  of  this  solid  unit,  with  oxygen  consumption 
on  the  cathode  side,  hydrogen  consumption  and  water  pro¬ 
duction  on  the  anode  side.  Heats  of  reactions  are  computed 
and  corresponding  heat  sources  are  included  in  the  energy 
equation.  Electric  power  is  removed  from  the  system. 

3.2.  Equations 


Every  equation  considers  the  three  dimensions  of  space 
and  the  time  dimension  however  in  the  present  study  the 
steady- state  condition  is  supposed. 

The  mass  conservation  equation  includes  transport  phe¬ 
nomena,  diffusion  and  a  reactive  production  term. 

—  (pYi)  +  V  (pvYi)  =  -VJi  +  n  (1) 

ot 

The  diffusion  term  is  calculated  from: 

N-l  y p 

Ji  =  -J2pDijVYj-DT4—  (2) 

j=  1 

Diffusion  coefficients  are  calculated  from  the  Maxwell-Stephan 
equation  [5]: 


xoq 

Pi 


(V 


1 


j=hjfi 


XiXj 


In  the  momentum  conservation  equation,  gravitational  body 
force  was  neglected  (gases,  forced  convection,  small  di¬ 
mensions).  Pressure  and  stress  tensor  work  terms  are  cal¬ 
culated  with  an  additional  source  term  F,  for  a  simulated 


To  model  the  gas  diffusion  layers,  we  assumed  a  porous 
medium  between  the  cell  and  the  interconnects.  The  porous 
media  is  then  modeled  as  a  sink  in  the  momentum  equation: 


3  3 

E  D'duiivj  +  y]  Cij  ^  Pvma.gV j 
7=1  7=1 


The  energy  equation  is  computed  from  a  porosity  weighted 
composite  of  energy  equations  of  the  solid  and  fluid  phase: 


—  (YPfEf  +  (1  -  y)psEs)  +  V (v(pfEf  +  p)) 
ot 


keffVT 


+  (tv) 


3.3.  Reaction  rate 


In  the  present  case,  the  fuel  is  hydrogen  and  the  oxidant 
is  air  so  the  overall  reaction  is: 

H2  +  3O2  H2O  (8) 

The  reaction  rate  rn2o  depends  on  current  density  j: 


m2o  = 


J 


IF 


(9) 


A  simple  semi-empirical  reaction  rate  is  used: 
f^cell  =  f^Nernst  (^ohm  T  Rpo\)j 


-AG0  RT 

f^Ncrnst  —  :m  b  t  m 


IF 


IF 


p1/2  pu 

r()2  rH2 

^h2o 


(10) 


(11) 


The  polarization  resistance  can  be  calculated  from  measure¬ 
ments.  On  small  cells  tests,  fuel  is  in  excess  so  that  partial 
pressures  in  Eq.  (11)  could  be  considered  as  constant.  A  re¬ 
sistance  function,  depending  on  temperature  and  local  cur¬ 
rent  is  determined.  We  fitted  the  resistance  with  a  temper¬ 
ature  and  current  dependent  law.  The  numerical  values  are 
available  in  Table  1. 

flpol  =  A(T)jmr>  (12) 

A  =  a\  x  T2  +  a2  x  T  +  <23 
B  =  b\  x  T2  +  Z?2  x  T  +  £>3 


Table  1 

Numerical  values 


a 

b 

1 

4.016E— 04 

— 3.283E— 06 

2 

— 8.509E— 01 

5.717E— 03 

3 

4.399E02 

— 2.792E00 
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Fig.  3.  equivalent  electric  scheme. 


A  constant  ohmic  resistance  (R0hm  =  0-23  fi/cm2)  is  added 
to  the  resistance  measured  on  small  cells,  to  take  into  ac¬ 
count  the  gas  diffusion  layer  and  the  interconnect  electrical 
resistance. 

With  the  ASE  concept,  the  electrolyte  can  be  very  thin, 
in  our  case  between  5  and  10  jutm. 

In  this  case,  the  electrolyte  can  be  doped  throughout  by 
impurities,  causing  electronic  conduction  increase  and  leak¬ 
age  current.  Electronic  conduction  in  one  direction  allows 
ionic  conduction  in  the  other  [6].  This  implies  a  modified 
open  circuit  voltage.  The  phenomenon  is  added  in  the  model. 
The  circuit  was  modified,  adding  an  electronic  conduction 
path  for  the  electrolyte  [7].  This  configuration  is  shown  in 
Fig.  3.  Leakage  currents  depend  on  the  ionic  and  electronic 
conductivity.  Those  values  were  estimated  [8]  from  experi¬ 
mental  OCV  values  at  several  temperatures,  Eq.  (10)  leads 
to  (13): 

f^cell  =  Uqcv  ~  (Rohm  +  Rpo\)j  (13) 


At  open  circuit,  ionic  current  is  equal  to  electronic  current: 


f^OCV  —  f^Nernst  ^ionicjleak 


(14a) 


Jleak  — 


Uqcv 

^electronic 


(14b) 


3.4.  Boundary  conditions 

The  cell  exchanges  heat  by  conduction  with  the  test  set-up 
(metal  flanges),  the  setup  is  placed  in  an  oven  and  exchanges 
heat  with  the  environment.  This  exchange  is  driven  by  two 
phenomena:  natural  convection  within  the  gas  of  the  oven 
and  radiation  between  inner  surfaces.  In  the  laboratory  test 
setup,  unreacted  fuel  is  burned  right  outside  the  cell  with 
oxygen  from  the  environment.  This  post-combustion  zone 
creates  a  higher  temperature  heat  exchange  with  the  corre¬ 
sponding  side.  Calculating  combustion  rates  implies  a  fine 
meshed  zone.  In  the  present  model,  we  simplified  the  phe¬ 
nomenon  by  fixing  a  combustion  temperature.  A  more  de¬ 
tailed  model  with  a  meshed  post-combustion  zone  is  in  elab¬ 
oration  at  this  time  and  will  be  presented  in  future.  As  the 
meshed  volume  stops  at  the  outlet  of  the  stack,  we  could  not 


calculate  species  back-diffusion  through  those  boundaries. 
We  will  notice  consequences  of  this  assumption  further  on. 

3.5.  Numerical  aspects 

Each  of  the  previous  equations,  including  the  reaction 
rates,  are  solved  at  each  iteration.  The  mesh  was  created  with 
GAMBIT  preprocessing  software.  The  model  includes  inter¬ 
connects,  gas  diffusion  layers  and  the  cell.  The  entire  meshed 
volume  is  divided  in  about  100,000  cells  (see  Fig.  4),  with 
variable  sizes  depending  on  location.  Because  of  a  symme¬ 
try  plane  in  the  geometry,  only  half  of  the  volume  is  meshed. 
Mesh  dimensions  follow  real  experimental  sizes  except  for 
interconnects  thicknesses  (1  mm)  in  the  model  where  exper¬ 
imental  interconnects  are  0.75  mm  thick.  The  model  yields 
results  for  a  set  of  boundary  conditions  in  about  150  itera¬ 
tions  and  10  min. 

3.6.  Validation 

As  functions  for  the  reaction  rates,  species  consumption, 
heat  sources  and  sinks  were  defined,  we  verified  the  mass 
and  energy  balances. 

The  mass  balance  between  inlet  and  outlet  mass  fluxes 
amounted  to: 

y>m-E  Mout  =  emass.  £mass  <  1.5  x  10-12  kg/s 

(15) 

The  energy  balance  between  heat  fluxes  through  outer  solid 
surfaces,  energy  fluxes  through  inlet  and  outlets  by  transport 
phenomena,  heats  of  reaction  and  electric  power  amounted 
to: 

MoutHout  +  E  MinHin  +  C/ceii/  =  £energy, 

^energy  <  2.8  X  10-3  W 

I  =  fsjdS  (16) 

H  =  Especies  Af  A?  +  Si  Cp  d  T 

We  can  conclude  that  the  sub-routines  are  consistent  with 
conservation  equations  and  that  the  model  is  numerically 
verified. 


Fig.  4.  In-plane  and  vertical  cut  of  the  mesh. 
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4.  Results  and  discussion 

Characteristic  I  —  V  curves  have  been  simulated  to  com¬ 
pare  with  experimental  values  (Figs.  5  and  6).  Besides  the 


general  trend  of  the  curves,  we  noticed  a  shift  of  the  ex¬ 
perimental  curve,  of  about  —20  mV;  the  experimental  open 
circuit  voltage  is  still  lower  than  the  calculated  one.  We  ex¬ 
plain  this  phenomenon  by  omitting  in  the  model  to  calculate 


Current  density  (A/cm2) 

Fig.  5.  Cell  potential  depending  on  current  density,  H2  flux  360  ml/ min,  airflux21/min,  50  cm2  active  area. 
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Fig.  6.  Cell  potential  depending  on  current  density;  200 ml/min  H2,  A  =  2.5. 
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Fig.  7.  Current  density  on  the  cell  (A/m2)  @0.7  V  and  28.8  A;  @40% 
fuel  utilization. 
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Fig.  9.  Pressure  field  in  anode  side. 


back-diffusion  through  outlets,  which  modify  partial  pres¬ 
sure  near  the  outlet  and  thus  the  OCV.  This  diffusion  could 
lead  to  the  observed  loss. 

Comparing  small  cell  results  to  repeat  element  results, 
a  strong  difference  is  apparent.  Average  current  densities 
on  large  repeat  elements  are  only  half  of  those  on  small 
cells,  at  the  same  potential.  For  a  large  part,  this  is  as¬ 
signed  to  the  geometrical  effect.  As  the  partial  pressure 
of  hydrogen  decreases,  Nernst  potential  decreases.  Fig.  7 
shows  current  density  distribution  for  the  repeat  element 
fed  with  360ml/min  hydrogen  and  an  air  excess  ratio  of 
2.5.  Near  the  fuel  entrance  (Fig.  7),  the  large  cell  pro¬ 
duces  locally  around  1  A  cm-2,  which  can  be  compared  to 
small  cells  results.  Large  areas  close  to  the  walls  encounter 
low  hydrogen  partial  pressure  due  to  lower  local  flow  rate 
(Fig.  8). 

Adequate  geometrical  design  is  therefore  paramount.  To 
increase  average  current  density,  effort  has  to  be  put  in  cell 
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Fig.  8.  Hydrogen  molar  fraction  @0.7  V  and  28.8  A;  @60%  fuel  utiliza¬ 
tion. 


and  stack  shaping.  In  this  case,  a  strong  limitation  of  perfor¬ 
mance  arise  from  geometry:  this  limitation  is  not  explained 
by  electrochemical  performance  but  by  a  non  homogeneous 
supply  of  reactants  on  the  surface. 

The  fine  mesh  of  the  model  allows  to  compute  the  pressure 
field  so  as  to  predict  pressure  drop  for  complex  manifold,  de¬ 
pending  on  temperature  gas  velocity  and  composition.  Fig.  9 
shows  an  example  of  the  pressure  results  obtained. 


5.  Conclusion 

A  modeling  tool  has  been  developed  and  validated  to 
compute  complex  flow  fields  in  SOFC  repeat  elements,  with 
detailed  geometry.  Good  accuracy  has  been  demonstrated 
between  simulations  and  experimental  results. 

The  methodology  employed  is  as  follows. 

(1)  Electrochemical  characterisation  of  ASE,  defining  the 
performance  within  a  range  of  temperature  and  current. 
This  experimental  work  is  done  on  small  cells. 

(2)  Modelling  the  real  geometry  of  the  large  cell  design  to 
predict  power  and  behavior. 

(3)  Optimize  and  modify  the  geometry  in  order  to  increase 
performance,  and  avoid  poor  concentration  regions. 

An  explanation  of  the  lower  power  density  obtained  on  a 
repeat  element  compared  to  a  small  cell  was  given.  These 
results  have  shown  that  cell  performance  can  be  improved 
by  better  flow  design.  The  need  for  optimization  in  the 
geometrical  domain  arises,  with  multiple  criteria  such  as 
power  density,  temperature  gradients  or  chemical  stability. 
In  further  developments,  the  post-combustion  zone  will  be 
simulated  to  study  the  interaction  with  the  cell  such  as  heat 
exchange  or  species  back-diffusion  through  the  outlet.  No 
further  development  of  the  model  is  needed  for  studying 
transient  behavior  as  the  software  is  already  designed  for 
it.  As  the  type  of  fuel  is  a  key  issue  for  fuel  cells,  internal 
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reforming  aspects  of  various  fuels  will  be  studied  adding 
subroutines  for  chemical  behavior. 
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